%% term_F_Psi
function vector_b = term_F_Psi(mesh, fespace, func_f)
Nb = fespace.Nb;
Nlb = fespace.Nlb;
rows = zeros(mesh.Ne*Nlb,1);
val = rows;
flag = 1;
for n = 1:mesh.Ne
    V = mesh.node(mesh.elem(n,:),:);
    [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
    phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace.basis_type, "x");
    ff = func_f(int_points_loc(:,1), int_points_loc(:,2));
    for i = 1:Nlb
        f = ff .* phi(:,i);
        r = int_weights_loc * f;
        rows(flag) = fespace.elem(n,i);
        val(flag) = r;
        flag = flag + 1;
    end
end
vector_b = sparse(rows, 1, val, Nb, 1);
end